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ABSTRACT 


We  perform  the  rounding  error  analysis  of  the  conjugate  gradient  algo- 

* ^ 

rithms  for  the  solution  of  a large  system  of  linear  equations  Ax  * b where 
A is  an  hermitian  and  positive  definite  matrix.  We  propose  a new  class  of 
conjugate  gradient  algorithms  and  prove  that  in  the  spectral  norm  the  rela- 
tive error  of  the  computed  sequence  (in  floating  point  arithmetic) 

depends  at  worst  on  ,£k  / where  £ is  the  relative  computer  precision  and  k 

i~  ?<- 4 «■ 


is  the  condition  number  of  A.  We  show  that  the  residual  vectors  r 


k * Avb  . 

- iu  * \C 


are  at  worst  of  order  £x||  A ||  ||  xk)|  . We  point  out  that  with  iterative  refine- 

2 

ment  these  algorithms  are  numerically  stable.  If  £n  is  at  most  of  order 
unity,  then  they  are  also  well-behaved.  . . 7 
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1.1 

1.  INTRODUCTION 

We  study  conjugate  gradient  algorithms  (for  brevity,  eg  algorithms) 

— * 

for  the  solution  of  a large  linear  system  Ax  ■ b where  A is  an  nxn  hermit- 
ian  and  positive  definite  matrix.  By  a eg  algorithm  we  mean  an  implementa- 
tion of  the  eg  iteration  in  floating  point  arithmetic.  We  are  primarily 
interested  in  the  round-off  error  analysis  of  these  algorithms. 

It  is  well-known  that  the  eg  iteration  enjoys  optimal  complexity  in  a 
sense  to  be  made  precise  in  Section  2.  In  exact  arithmetic  it  generates  a 
sequence  of  orthogonal  residual  vectors  r^  * Ax^-b  anc*  s°luti°n  a m A 
is  obtained  after  at  most  n steps.  See,  among  others,  Hestenes  and  Stiefel 

[52],  Stiefel  [58]  and  Engeli,  Ginsburg,  Rutishauser  and  Stiefel  [59]. 

■ 

Many  of  these  theoretical  properties  do  not  hold  in  the  presence  of  rounding 

errors.  It  is  no  longer  true  that  the  computed  residual  vectors  are  ortho- 

- 

gonal  (or  even  nearly  orthogonal)  and  that  the  nth  computed  vector  xn  is  a 
reasonable  approximation  to  a. 

The  aim  of  this  paper  is  to  understand  the  behavior  of  some  eg  algorithms 

in  the  presence  of  rounding  errors.  We  are  primarily  interested  in  studying 

“1  . 

how  the  matrix  condition  number  x =■  ||  A j j ||  A j|,  where  ||  A J | denotes  the  spec- 
tral norm  of  A,  influences  the  relative  error  of  the  computed  sequence  {x^}. 

We  know  that  direct  algorithms  of  practical  interest  as  well  as  many 
iterative  algorithms  with  iterative  refinement  are  well-behaved,  i.e.,  they 
assure  the  computation  of  an  approximation  y in  floating  point  arithmetic  fl 
I . such  that  y is  the  exact  solution  of  a slightly  perturbed  system,  i.e., 

I(A+6A)y  ■ b where  ||  6A  |j  is  of  order  £||  A||  and  £ is  the  relative  computer 

precision.  Equivalently,  the  residual  vector  r ■ Ay-b  has  the  norm  of  order 


1.2 


C||  A ||  ||  y ||  . When  it  cannot  be  established  that  an  algorithm  is  well-behaved 
it  is  sometimes  possible  to  prove  a weaker  property,  namely  that  an  algorithm 
is  numerically  stable,  i.e.,  the  relative  error  of  a computed  y is  of  order 
£h.  See  Wilkinson  [63  and  65]  for  direct  algorithms,  and  Jankowski  and 
Wozniakowski  [77],  Wozniakowski  [77  and  78]  for  iterative  algorithms. 

To  help  the  reader  get  an  idea  of  what  results  can  be  expected  for  eg 
algorithms  we  report  numerical  tests.  We  tested  the  Concus,  Golub  and 
O'Leary  [76]  conjugate  gradient  algorithm.  We  performed  j,  j » n,  iterative 
steps  finding  the  best  possible  approximation  x^,-  k £ j,  among  all  computed 
vectors.  Next  we  computed  the  relative  error  of  and  its  residual  vector. 
Define  the  number  s such  that 
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1.3 


Note  that  the  condition  number  of  A in  the  A-norm  coincides  with  the  condi- 
tion number  of  A in  the  spectral  norm.  Since  ||  x^-al)  * ||  A ||  ||  A1^2 ||» 


(1.2)  yields 


II  V“H  - °(c»3/2ll  > 


This  explains  why  s * 3/2  might  be  expected  in  (1.1).  See  the  Appendix  where 
a detailed  discussion  of  numerical  tests  is  reported. 

We  have  not  succeeded  in  analyzing  classical  eg  algorithms  including 
that  proposed  by  Concus , Golub  and  O'Leary  [76].  In  this  paper  we  propose  a 
new  class  of  eg  algorithms  and  prove  that  for  these  algorithms  there  exists  a 
computed  vector  such  that 


(1.3)  II  A1//2(^-aj  ||  s CCh||  a1/2||  II  ^11 


where  C is  a constant  of  order  at  most  n.  We  shall  denote  this  class  of  algo- 
rithms by  $.  Note  that  h occurs  linearly  in  (1.3).  In  general,  we  cannot  say 
that  (1.3)  means  numerical  stability  of  the  eg  algorithms  in  its  "own"  norm 


since  we  have  A 


x^||  instead  of  ||  A^2XjJ|  . However,  if  ||  A^2  | 


is  of  order  ||  A XjJ|  then  these  eg  algorithms  are  numerically  stable  in  the 


A-  no  t iii  • 


For  the  residual  vectors  we  are  only  able  to  prove  that 


(1.4)  ||  ?k  ||  SCC*  ||  A||  ||  ^|| 


We  tested  one  algorithm  cp  from  §.  For  most  cases  cp  looked  like  a well-behaved 
algorithm,  i.e.,  ||  r^||  was  of  order  C||  A||  ||  x^||  . However,  for  a few  cases, 

||  r^||  was  of  order  £*||  A 1 1 ||  x^  ||  . This  proves  that  (1.4)  is  sharp  and  some 
eg  algorithms  from  * are  not  well-behaved. 


1.4 


Many  iterative  algorithms  have  this  property,  i.e.,  they  are  numerically 

stable  but  not  well-behaved.  Examples  include  the  Chebyshev,  SOR,  Richardson 

and  Jacobi  iterative  algorithms.  See  Wozniakowski  [77  and  78].  However,  it 

was  shown  by  Jankowski  and  Woz'niakowski  [77]  that  any  algorithm  (direct  or 

iterative)  which  computes  an  approximation  y such  that  ||  y-o||  £ q ||  a||  with 

q < 1 followed  by  iterative  refinement  in  single  precision  becomes  numerically 
2 

stable  and  if  ■ £*  is  of  order  unity  then  it  is  also  well-behaved.  Hence,  any 

eg  algorithm  from  5 with  iterative  refinement  is  numerically  stable  in  the 

3/2 

spectral  norm  whenever  is  bounded  away  from  unity  and  it  is  well-behaved 

2 

whenever  is  of  order  unity. 

We  summarize  the  contents  of  the  paper.  In  Section  2 we  briefly  state 
the  basic  theoretical  properties  of  the  eg  iteration  and  derive  a three  term 
recurrence  formula  for  the  vectors  [x^}  which  explains  the  connection  between 
the  eg  iteration  and  the  steepest  descent  iteration. 

Section  3 deals  with  the  rounding  error  analysis  of  a steepest  descent 
algorithm.  We  prove  that  inequality  (1.3)  holds  for  this  algorithm. 

Section  4 deals  with  the  round-off  error  analysis  of  a new  class  of  eg 
algorithms.  Based  on  the  results  of  Section  3 we  prove  (1.3)  and  (1.4)  which 
are  the  main  results  of  this  paper. 

In  the  final  section  we  pose  a conjecture  on  the  speed  of  convergence  of 
sequences  computed  by  eg  algorithms. 


2.1 


2.  GRADIENT  AND  CONJUGATE  GRADIENT  ITERATIONS 

In  this  section  we  briefly  derive  some  basic  properties  of  the  gradient 
and  conjugate  gradient  iterations.  We  consider  the  solution  of  a large  linear 
system 

(2.1)  Ax  * b 


where  A = A > 0 is  an  nxn  hermitian  and  positive  definite  matrix  and  b is 
a nxl  given  vector.  Suppose  that  the  information  about  the  matrix  A is  given 
by  a procedure  which  computes  y = Ax  for  a given  x.  For  large  systems  A is 
usually  sparse  which  permits  the  evaluation  of  y in  time  and  storage  propor- 
tional to  n. 

We  solve  (2.1)  iteratively  by  constructing  a sequence  [x^}  converging 
— - 1-*  * 

to  the  solution  a * A b.  Let  B = B > 0 be  a matrix  which  commutes  with  A, 

BA  * AB.  For  instance  one  can  set  B = AP  for  a real  p.  Let  ||  x||g  “ V(Bx,x)  ■ 
* ||  B^2x||  where  ||  x||  =»  «/(x,x)  is  the  spectral  norm. 

We  recall  the  definition  of  the  gradient  iteration  which  constructs  the 
sequence  [x^l  as  follows.  Let  x^  be  a given  initial  approximation  and 


(2-2)  Vi  = *k  • Vk’  rk  = Avb> 

where  c,  is  chosen  in  such  a way  that  the  error  e,  , 
k k+1 

i-e*>  II  xk+r®llB  = inf  II  V^klt'  This  yields 
c 

<*k,B(*k~“)) 


x^+^-o||g  is  minimized, 


(2.3)  c. 


<rk>Brk) 


Note  that  (B(x^+1-aj , r^)  = 0.  The  coefficient  c^  is  computable  only  for  certain 


2.2 


matrices  B.  Suppose  that  B 3 AP  where  p is  an  integer.  Then 

= (rk’AP  comPutable-  For  B = I,  we  cannot,  in  general, 

compute  the  numerator  of  (2.3),  (r^,x^-or)  . However  if  one  considers  the  system 
Mx  = g with  a nonsingular  M which  is  nonhermitian  or  nonpositive  definite  and 
if  one  agrees  to  multiply  this  system  by  M then  A = M^M,  b = M*g  and 

Ck  = (Mxk"S> M^ic'S)/ (rk>rk)  iS  comPutable« 

It  is  well-known  that  {x^}  converges  to  a and 

k+1 


(2.4)  «ktl  - v4^-ck(Jk,Brk)  s 


where  e^_  = x^-  a,  e^  - ||  e^||g  and  a - 1 1 A (|  ||  A ^ (|  is  the  condition  number  of 
matrix  A.  (Note  that  ||  A^  = ||  A }]  and  ||  A-1^  =>  ||  A~ 1 J|  .) 

Recall  that  for  B = A,  (2.2)  and  (2.3)  is  called  the  steepest  descent 
iteration.  It  has,  in  general,  very  slow  convergence  and  therefore  is  not 
recommended  in  numerical  practice.  The  conjugate  gradient  iteration  is  much 
more  efficient.  The  following  derivation  of  the  eg  iteration  focuses  on  its 
complexity  optimality. 

Consider  a class  of  iterations  for  which  the  error  formula  satisfies  the 
relation 


(2.5)  xk-o<  = Wk(A)  (xQ-of) 

where  Wk  is  a polynomial  of  degree  at  most  k and  Wk(0)  = 1.  A natural  complex- 
ity question  is  how  to  choose  the  polynomials  Wk<  Since  we  want  to  minimize 
the  computational  complexity  (cost)  we  seek  wk  so  that  the  error  efc  « || 
is  minimized.  This  means  that  the  polynomials  V?k  are  the  solution  of  the  fol- 
lowing problem 
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k. 


1 


m 

(2.8)  (f,g)  = Kl  |Cj|2  f3j Xjf (XJ)iTxT) 

J-l 


where  f and  g are  functions  defined  on  the  interval  The  polynomials 

Wk>  Wk(0)  = 1>  which  minimize  (2.6)  are  the  orthogonal  polynomials  with  respect 
to  the  inner  product  (2.8),  i.e., 
m 

(2.9)  (Wk,W.)  - l |Cj|2  Q.\.  Wk(X.)Wi(X.)  =•  0 

for  k / i.  From  the  orthogonality  of  Wk  it  follows  that  they  satisfy  a three- 
term  recurrence  formula.  We  choose  a different  form  of  the  three- term  recur- 
rence formula  than  usual  in  order  to  emphasize  the  connection  between  the  eg 
iteration  and  the  gradient  one.  This  form  is  defined  as  follows: 


2.4 


W„(X)  = 1 


(2.10)  W1(X)  - 1 - cQX 

Vi(W  ■ Vx>  • ckwk<w  - ”kl"k-i(X)-"k<x>+  'i'l™1'  k 2 1 


where 


Cw  ,w  ) 

(2.H)  Ck=(AWk,Wk)’ 


(2.12)  u0=0,  uk= 


(wk_ck^k*  \(wk-  rV+ckV 


k s i. 


(\-  fwk+\M,k’\(wk-  rwk)+ckMk} 

From  this  we  get  the  three- term  recurrence  formula  for  the  sequence 


z.  = x,  - c.  r,  , 
(2.i3)  _k  k_  kk: 

Vi  = zk  - Vk- 


-b. 


rk  = Axk' 


yk  = Vi  - V 


From  (2.11),  (2.12)  and  (2.9)  we  get 


(rk*  B(va)) 


(2.14) 


(V  ^ 


0.  ub 


(yk>  B(vcy)) 

(v  BV 


k a 1. 


The  conjugate  gradient  iteration  (2.13)  consists  of  computing  z,  and  x 

k k+1 

The  vector  zk  is  obtained  by  one  step  of  the  gradient  iteration  (2.2)  and 

is  the  best  approximation  of  y along  the  line  rk-  The  vector  x^^  is  the 

best  approximation  of  y along  the  line  yk,  ||  x^+^_°i!g  “ inf  II  zk_IrU^k^' 

u 

From  the  orthogonality  of  Wk  it  follows  that 


(2.15) 


(B(xk-o),  rj)  = 0 
(B(xk-a),  3^-x  ) = 0 


for  j * k, 
for  j s k 


2.5 


This  yields  = ck(rk,  B^-a) )/ (yfc,  Byk)  which  can  be  computed  for  B = I when- 
* 

ever  A =MM,. 

The  conjugate  gradient  method  (2.13)  converges  in  exactly  m steps,  i.e., 

(2.16)  = ni  for  k s m. 

From  (2.6)  one  can  estimate  the  speed  of  convergence  for  initial  approxima- 
tions xk>  k £ m.  Setting  P(\)  = Tk(f (X) )/Tk(f (0) ) 'in  (2.6)  where  Tk  is  the 

kth  Chebyshev  polynomial  of  degree  k and  f(X)  = (X +X.. -2  X)/ (X  - X ) we  get 

ml  ml 

(2.17)  ||  £ ||  P(A)  (xQ-a)  ^ * 2^‘^)  ||  xQ-o|^ 


where  a = X^/X^.  Compare  with  (2.4).  For  the  spectral  norm  we  have 

<2.ls)  iiv*Hs^(^r)kllv3il 


3.1 


i 


3.  ROUND-OFF  ERROR  ANALYSIS  OF  GRADIENT  ALGORITHMS 


We  shall  show  that  the  round-off  error  analysis  of  eg  algorithms  belonging 
to  $ can  be  primarily  based  on  the  round-off  error  analysis  of  the  gradient 
algorithms  to  be  studied  in  this  section.  Therefore  in  this  section  we  analyze 
gradient  algorithms  in  the  presence  of  rounding  errors.  We  focus  our  atten- 
tion on  a steepest  descent  algorithm  (B  =*  A)  and  mention  the  corresponding 

2 

results  for  the  gradient  algorithms  with  B = I or  B “ A . 

We  consider  a steepest  descent  algorithm  in  floating  point  binary  arith- 
metic, fl,  with  the  relative  computer  precision  £ m 2 t where  t is  the  number 
of  mantissa  bits.  To  simplify  further  estimates  we  shall  use  the  relation 
which  is  defined  as  follows.  Let  f and  h be  two  scalar  functions  defined 
on  [0,Cq].  By  f(£)  * h(£)  we  mean  that  there  exists  a constant  c such  that 
f(0  = h(Q(l  + e(O)  where  |e(Ol  4 cC  for  0 s C s £Q.  By  f(£)  f h(£)  we 

mean  f (£)  s h(£)  or  f(£)  ^ h(£) . The  relation  enables  us  to  ignore  the 
2 

terms  of  order  £ in  the  presence  of  the  term  of  order  £. 


Let  and  r^  denote  the  vectors  computed  in  fl  by  an  algorithm.  We 


as  sume  tha  t 


(3.1)  rk  ■ f 1 (Ax^-b) . 

Let  r^  = Ax^-b  denote  the  exact  value  of  the  residual  vector  and  let 

(3.2)  ek  * A (lyor),  eR  = ||  ek||  . 


We  analyze  (3.1).  Assume  that  the  algorithm  for  evaluation  of  rk  satisfies 
the  relation 


3.2 


(3.3)  rk  - (I+D^)  ( (A+E^)  x^-b)  - r*  + 6?k 


where  D*  is  a diagonal  matrix  such  that  ||  Dk||  £ Q and  ||  Ek||  * £||  A ||  C1  where 
the  constant  depends  only  on  the  size  of  the  problem,  “ C^(n).  Hence 


.1- 


1,-*  1- 


(3.4) 


ck  ■ Ek\  + Dk(tk + EkV 


II  «?kll  fell  All  ||^l|C1  + dl^|| 


-a* 


To  assure  that  r,  is  a reasonable  approximation  to  r,  we  assume  that  ||  r ||  > £||  a| 


II  XjJI  C Note  that  the  opposite  inequality  ||  r^ ||  < C||  a||  ||  XjJI  means  that 

is  the  exact  solution  of  the  system  (A-SAjx^^  * b where  5A  ■ (rk*  5rk)x£/ 1|  XjJp 
and  ||  «A||  ^ C ||  A||  2C^.  This  means  that  the  algorithm  is  well-behaved  and  the 
iteration  should  be  terminated.  Therefore  we  shall  analyze  the  following  algo- 
rithm of  the  steepest  descent  iteration. 


Algorithm  3.1 


Let  Xq  be  given  and  let  k ■ 0. 


(i)  compute  rk  = fl(Axk~b), 

(ii)  if  ||  rk||  ^ C IJ  A ||  ||  xk||  then  formally  set  Xj^  • xk>  Vi  a k+1. 
STOP. 

(iii)  if  II  rkll  > C ||  A ||  ||  xk||  Cl  then  compute 


(3.5)  ck  = fl((rk,rk)/(rk,  Ark)), 


*k+l  * fl(Vckrk)’ 

(3.6)  k :»  k+1 
GO  TO  (i) 


We  do  not  define  a termination  criterion  for  Algorithm  3.1  unless  it  is 
well-behaved.  We  want  to  verify  its  numerical  stability  by  Tim H^- or||  and 


V 


has  to  be  defined  for  all  k.  Therefore  we  formally  set  x^  = x^,  Vi  2 k+1 
in  (ii) . 

Recall  that  the  computed  inner  product  (a,b)  in  fl  satisfies  the  relation 
(3.7)  fl(a,b))  = ((I+D)a,b) 

where  D is  a diagonal  matrix  such  that  ||  d||  £ £c2.  The  constant  C2  depends  on 
a particular  algorithm  used  for  the  summation  of  n numbers.  For  the  standard 
algorithm  = n whereas  when  the  Holier  algorithm  is  employed  C2  y 3.  See 
Wilkinson  [63]  for  the  first  and  Kielbasinski  [73]  for  the  second  result. 

We  are  ready  to  prove 


Lemma  3 . 1 

Suppose  that  j|  r^||  > C||  A 1 1 ||  x^j|  C^,  Vk.  Then  the  sequence  computed 

by  Algorithm  3.1  satisfies  the  following  error  formula: 


<3-8>  Vi  f'Vck 


Cif 


III  A1-'2, 


* 


,3/2 


A'  'll  II  \W  5C: 

A||  ek(C1+2C2+8)} 


where  r, 


Axk'^J  Ck 


(7k’7k)/(h’AK} 


Proof 

We  analyze  the  computation  of  c^.  Due  to  the  assumption  that 
||  tjJI  > s||  A ||  ||  XjJI  C^,  Vk,  (3.3)  and  (3.4)  imply  that  r^  * 0 and  c^  is 
well  defined.  We  have 


(a+Dk)r  ,r  ) j 

(3.9)  c - 3 2--  (1+V 

((I+Dk)(A+Ek)rk,rk) 


(ru>rJ  2 

(l+ek) 

(rk'Ark> 


i 2 1 

where  |j  DR  ||  * £c,  for  i = 2 and  3,  ||  Ek||  * ;||  a||  ek  is  the  relative  error 


of  division,  |e^|  ^ Q,  and 

✓ ,_2-  - 


2 m ( , i°kVvY  „ (CEH(WEfc)  ]vv 
k v TwTA 


(3.10)  , / 

l\l  f C(C2  + 1 + ||  A||  (C1+C2) 


w 

(rk,Ark) 


From  (3.3)  we  get 


(rk  ‘V-  * 3 

— — " ck(1+ek>> 

(rk,Ark) 


3|  2||6rk||  2||AJ 


(3.11) 


f 2C(||  A||  ||^||C1+  ||  rk||  ) -(-i 


A 1/2—* 
A ' T 

A rk 


f4C| |A1/2||  (||  A||  ||  ^||  Cx 


ai/2?:i 


From  (3.9)  and  (3.11)  we  have 


Ck  = Ck(1+6ck)‘ 

(3.12)  |6ck|  f l\l  + l\l  f C[C2  + l + c*||a||  (C^)  + 


4 II  A1/ 2 II  (II  All  |f^||c1+  H^||)/||A 


We  now  analyze  (3.6).  We  have 


(3.!3)  xk+1  “ d+Dk)  (\-(I+Dk)Ckrk)  “ xk  “ Ckrk  + ^+1 


where  ||  Dk||  s C for  i ■ 4 and  5, 


3.5 


(3.14)  5^  = 


* -*  * -kiv  4 5 -*  2 

- ck5rk  - ck*5ck*rk  - Dkckrk  + ^(x^-c^)  + 0(£  ) 


From  (3.4),  (3.12)  and  (3.14)  we  get 


* 1/2—*  1/2  - 

6k+l  = ek  - CkA  rk  + A 5xk+l> 


(3.15) 


Al//2*xk+1ll  f C||  Al/2!1  ||xj|  + Cck[(|A3/2jj  ||xj|  5C1 

+ II  A II  ek(C2+8)  ] + 

+ t(<>2||*||  ll*1/2^li  (ct+c2) . 


Note  that 


<11  Al/2?kH  = <Al/2^»«k)/H  A*  2\il  s V 


From  this  and  (2.4)  we  get  (3.8)  which  completes  the  proof  of  Lemma  3.1.  ■ 

Lemma  3.1  shows  how  the  error  ek+^  depends  on  the  theoretical  and  rounding 
errors.  It  is  interesting  to  notice  that  the  bound  on  the  rounding  error 

•k 

increases  with  ck  whereas  the  bound  on  the  theoretical  error  decreases  with 
* 

increasing  ck« 

We  want  to  find  the  limiting  properties  of  the  sequence  £&k } which  satis- 
fies (3.8).  To  achieve  this  we  use  the  following  lemma. 


Lemma  3 . 2 


2 *n-*(|2  * * 

ek+i  * vvckH  rkir  + \ + ck\  + ckekd 


for  given  nonnegative  sequences  Uk},  [bk}  and  a constant  d such  that 

2d  ||  a"1  II  < 1.  Then 


(3.16)  lim  ek  s 3x(lim  t^/  ||  A ||  + lim  ak)/(l  - 2d»t/ 1|  A ||  ). 
k k k 


4 


3.6 


Proof 

Let  e be  any  positive  number.  Choose  k^  such  that  a^-e  £ a m lim  a^, 

* k 

b -e  £ b ■ lim  b^  for  k ^ kQ.  Then  ek+^  £ f(c^)  w^ere 

k k 


f(c)  - vVcll;klf  + (a+e)  + c(bfe)  + cekd. 


for  c 6 [ ||  A ||”\  ||  A 1||  ].  Consider  two  cases. 

Case  (i).  Assume  that  e^  s 2 ||  A 1 1|  (b+e)/  (1  - 2 ||  A 1 1|  d) , Vk  ^ kQ.  We  show 
that  f is  decreasing  for  c S 0.  Indeed, 


- II  r* 


f'  (c) 


If 


- ii  ?,* 


2"<VC  II  ?£lf 
ii?;  if 


+ b + e + dek  £ f'  (0) 


If 


. + b + e +■  de,  . 

2e,  k 


kk 11  / 1 

Since  de, s £ I d r-  1 e. 


'k  2e, 


2 A 


-In  k 


1-  (l-2d  ||  A~ 1 1|  ) e,  £ - (b+e),  f' (c)  *0 

2||  A'1))  k 


we  get  f'  (c)  £ 0.  Thus 
Since  ||  r*|p/||  A ||  £ ^I-h”1  ek  we  get 


f(c)  £ f (||  Alf1)  - v4k-||  r*|f/||  A||fa+e+(b+e)/||  A|^d/||  A j(  . 


-1 


(3.17)  e . £ «/l-K  e + (a+e)  + (b+e)/||A||+  e d||A||  . 


Note  that  [e^}  is  a nonincreasing  sequence  and 


lim  e,  £ (a  + e + (b+-e)/ ||  A ||  )/ (1  - */l-H  1 - d/||  A ||  ) 


Since  1 - Ji-  1 * i/  (2k)  we  finally  get 


(3.18)  lim  £ P0  — 2*(a  + e + (b+e)/||  A j|  )/ (1  - 2d*/||  A ||  ) . 


3.7 


Case  (ii)  . Assume  there  exists  k,  such  that  e,  < 2 ||  A * II  (b+e)/ (1-2  ||  A ||  d)  £ 3 

i ^ 

We  inductively  prove  that 
(3.19)  ek  s 1.5  3qJ  k * k^ 

This  holds  for  k = k^.  Suppose  first  that  e^  satisfies  (3.19)  and  addition- 
ally e^  ^ 3g.  From  (3.17)  we  get 

6k+l  * vi-x."1  ek  + a + e + (b+s)/||  a||  + ekd/ ||  A ||  < \ * 1-5  3Q. 

If  e.  < then 
k 0 

ek+l  s ek  + a + e + (b+c)  ||  A'^l  + ||  A" 1 1|  ekd  * 

(1  + dx/  ||  A ||  ) + (a  + e + (b+e)/ ||  A ||  )h  - SQ(1  + d>./  j|  A jj 
+ .5(1  - 2dy./\\  A||  ))  - 1.5  3q. 

Hence,  in  all  cases  we  proved  that-  lim  e s 1.5  3 . Letting  e tend  to  _ero 

krC  U 

we  get  (3.16).  This  completes  the  proof.  ■ 

From  Lemmas  3.1  and  3.2  we  immediately  conclude  the  asymptotic  behavior 
of  the  sequence  (.x  } computed  by  Algorithm  3.1. 

Theorem  3.1 

If  3 « 2 C H(C1  + 2C2  + 3)  < 1 then  Algorithm  3.1  computes  the  sequence 
such  that 

— _ 1/2  _ _ 3(5C.  + 1)  ,/» 

(3.20)  lim  ||  A ' ‘(y3)  II  f & TT?-  II  A 


3.8 


4 


Proof 


Suppose  first  that  there  exists  kn  such  that  ||  r,  ||  £ £||  a||  ||  x,  ||  C, . 

U Kq  ICq  1 

Then  Algorithm  3.1  yields  x.  = x for  i s k . From  (3.3)  and  (3.4)  we  get 

V kQ  0 


lim  II  A1//2(^-o)  ||  = ||  A1//2(^o-a)  ||  ^ Ha'1/2)! 

Ch1^22c1||  A1//2||Tto  ||  xj| 


which  obviously  proves  (3.20). 

Hence  we  can  now  assume  that  ||  r.J|  > C||  A||  * ||  ^11  c^»  Vk.  Applying 
Lemma  3.2  with  ak  = Cll  A^2  ||  • II  *iJI  » bk  - C||  A3'  2 ||  ||  XjJI  5C1  and 
d = C||  A||  (C^+2C2+8)  we  get  (3.20)  from  Lemma  3.1. 

Theorem  3.1  states  that  if  k is  large  and  ||  x^U  — lim  ||  x^H  then  the 
— — k 

computed  x^  approximates  a with  the  error 


(3.21)  ||  Al/2(^-a)||  * Gk||  Al/2||  ||  C 

where  C = 3(5C^+1)  + 0(C) • Note  that  (3.21)  does  not  imply  the  numerical 

l/o  1 / 0— » 

stability  of  Algorithm  3.1  since  we  have  ||  A ' ||  ||  x^U  instead  of  ||  A XjJ|  . 

From  (3.21)  we  get 


(3.22) 


llv«llA  ( II  *1/2  II  II 


C. 


This  means  that  the  relative  error  of  x^  in  the  A-norm  depends  at  worst  on 
Ck3/  2 . However  if  ||  A*7 2 ||  x^ ||  a:  ||  A^"7  2x^ ||  then  Algorithm  3.1  is  numerically 
stable  in  the  A-norm. 

We  pass  to  results  for  the  spectral  norm.  From  (3.21)  we  have 


r 


3.9 


II  V«ll  II  a*1//2|I  II  Al/2(va)  II  3/2 

(3.23)  £ * ^ C. 


11  \ I! 


However  if 


(3.24)  ||  Al/2(^-;)||3=  ||  Al/2 


then  we  get  numerical  stability  in  the  spectral  norm.  Note  that  (3.24)  will 

n 

often  hold.  For  instance,  let  x.-ct  = c f where  5-  are  the  eigenvectors 

K — J J 3 

j=1  ,/2  _ _ 

of  A.  Suppose  that  c.  = c (or  c . 2=  c)  for  all  j.  Then  ||  A ' (x,  - ?)  j|  “ 

J J K 

/ - \ 1/2 

II  A1  2|l  ^i/\nax)2  J and>  II  A1  "II  llxk"^li=  II  A1  2 ||  |c|./n  and  these 

two  quantities  differ  at  most  by  a factor  of  n and  (3.24)  holds. 


For  the  residual  vector  r^  * Ax^  - b we  get 


(3.25)  ||rk||s;>. 


//2iD%5V,NI  !i^l|cs:K||Ail  l|Jkl|c- 


Numerical  tests  confirm  that  the  residual  vectors  sometimes  depends  ou  £*. 

This  means  that  Algorithm  3.1  is  not  well-behaved.  However  if 

||  A^-b  ||  — ||  A^  ^ (x^-or)  ||  / 1|  A ^ ||  then  the  residual  vector  r^  depends  at 

, 1/2 

worst  on 

Numerical  stability  and/or  well-behaved  property  may  be  achieved  by  the 
use  of  iterative  refinement  even  if  the  residuals  are  computed  in  single  pre- 
cision. From  Theorems  3.1  and  4.3  in  Jankowski  and  Wozniakowski  [77]  it 

follows  that  Algorithm  3.1  with  iterative  refinement  in  single  precision  is 

3 2 2 

numerically  stable  whenever  ' ~C  < 1 and  it  is  well-behaved  whenever  J* 

is  at  most  of  order  unity. 


3.10 


We  summarize  the  properties  of  Algorithm  3.1. 

Corollary  3.1 

Algorithm  3.1  constructs  an  approximation  x^  such  that 

IU1/2<V»)||scMI *l/2ll  II ^ He, 

II  vS||  s c«3/2  II  ^ II  c 

II  Av^H  * £*11  AH  ll  *kll  c 

where  C * 3(5C^+1)  + 0(Q  . Furthermore  if  ||  A1^2  ||  j)  XjJj  — jj  a^2x^ |j  then  the 
algorithm  is  numerically  stable  in  the  A-norm  and  if  jj  A ; ” ||  x^-ajl  — jj  A 7 (x^-a)  || 
then  the  algorithm  is  numerically  stable  in  the  spectral  norm.  ■ 

Corollary  3.1  summarizes  the  numerical  properties  of  the  steepest  descent 
algorithm.  It  shows  that  the  algorithm  may  be  neither  well-behaved  nor  numer- 
ically stable.  However  the  algorithm  is  guaranteed  to  compute  an  approxima- 

3/2 

tlon  with  a relative  error  of  order  at  most  CH  . If  the  problem  is  not  too 

ill-conditioned  this  is  a satisfactory  result. 

We  end  this  section  by  a remark  concerning  the  gradient  algorithms  for 
2 

B » I or  B * A . They  differ  from  Algorithm  3.1  by  different  computations  of 
c^  in  (3.5).  Based  on  proof  techniques  similar  to  those  used  here  it  is  pos- 
sible to  show  that  there  exists  an  index  k such  that 

II  »1/2Sk-3>lliC‘‘ll  »l/:H  ll\llc. 

llv»H‘°‘ll b1/2 ||  || b' 1/2 1|  n;j|c, 

II  V“l!  s HI  All  ||  £J|  c 

for  a certain  C ■ C(n).  This  shows  that  the  best  estimates  are  obtained  in 
the  "natural"  norm  of  the  algorithm  (i.e.,  in  the  B-norm)  and  that  the  residual 
vectors  may  depend  on  £h  for  every  choice  of  B. 


r 


mam 
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4.  ROUND-OFF  ERROR  ANALYSIS  OF  A CLASS  OF  CONJUGATE  GRADIENT  ALGORITHMS 

We  deal  with  the  conjugate  gradient  iteration  for  B = A which  generates 
the  sequence  as  follows. 

Zk  = *k  ' Ck  V \ = Axk  ' S’ 

(4.D  _ _ _ 

xk+i  = zk  ' \ V yk  = Vi  • V 


where 


(4.2) 


(Vrk} 

■*  ss  

(vAV  * 


uo  “ °’  uk 


(yk,A(zk‘a)) 

(yk)Ayk} 


, k a 1. 


See  (2.13)  and  (2.14).  It  was  pointed  out  to  the  author  by  Wieladek  [77]  that 
(4.1)  has  an  interesting  local  property.  Namely,  no  matter  how  the  vectors 
and  y^  are  computed,  the  coefficient  u^  is  chosen  in  such  a way  that  the 
error  ||  x^+^“all^  *-s  minimized  along  the  line  y^.  Note  that  the  cost  of  one 
step  of  the  eg  iteration  depends  on  how  one  computes  the  residual  vectors  and 
the  coefficients  c^  and  u^.  The  number  of  matrix- vector  multiplications  needed 
to  perform  one  step  may  vary  from  one  to  four. 

We  define  a new  class  of  eg  algorithms  § by  the  following  properties.  We 
assume  that  any  algorithm  cp  from  the  class  ? computes  the  vector  z^  by  Algorithm 
3.1.  That  is 


(4.3)  rk  = fl(Axk-b) 


and  if  ||  rk||  > C||  a||  ||  xjlc  l then 


4.2 


(4.4)  ck  3 fl((rk,rk)/(rk,Ark)), 

(4.5)  \ = fl(^  - ck  ?k). 

Thus,  the  computation  of  z may  require  two  matrix- vector  multiplications. 

• There  are  many  different  ways  of  computing  the  coefficient  uk<  One  may 

use  theoretical  orthogonality  relations  (2.15)  as  well  as  the  direct  substi- 
tutions for  yk  and  zk  from  (4.1).  For  instance  it  follows  from  (2.15)  and  (4.1) 
that  in  theory  uk  > 0.  For  the  sake  of  generality  we  do  not  specify  an  algo- 
rithm for  the  computation  of  u,  . We  only  assume  that  an  algorithm  cp  computes 

1 

uk  such  that 

(4.6)  5k  - uk(l  + 6uk),  |6uk|  * 1 

where  uk  ■ (yk>A(zk-a) )/ (yk*Ayk)  for  the  computed  vectors  zk  and 
yk  * fl(xk_ ^ - zk) . Note  that  (4.6)  means  that  iT^  can  be  a very  crude  approxi- 
mation of  uk>  A particular  algorithm  for  which  (4.6)  holds  is  given  in 
Example  4.1.  Knowing  uk  we  finally  compute 

(4.?)  ^+1  " fK*,.  - Vk).  j 

Thus  the  class  $ contains  algorithms  which  differ  by  the  computation  of  u . 

iC 

We  are  ready  to  prove 

Lemna  4.1 

Let  <p  be  a eg  algorithm  defined  by  (4.3)  to  (4.7).  Then 
(4.8)  ||  \+l-\  * (1+2C)||  V«lk  + Cll  Al/2||  ||xk+1||/(l-Q.  ■ 


4 


4.4 


^ _ 3 ( 5 C,  +2 ) _ 

(4.11)  lim  ||  A ' (\~a)  ||  * C* — jTT^ lim  II  \j|  . ■ 

^ k 

Theorem  4.1  states  the  numerical  properties  of  the  eg  algorithms  from 
the  class  $.  Since  (4.11)  is  essentially  equivalent  to  (3.20),  the  discussion 
of  numerical  properties  of  the  steepest  gradient  algorithm  is  also  valid  for 
the  eg  algorithms  from  $.  In  particular  we  can  estimate  the  error  and 

the  residual  vector  rfe  in  the  spectral  norm,  as  in  (3.22)  to  (3.25).  This  is 
summarized  in  the  following  corollary. 

Corollary  4 . 1 

Any  eg  algorithm  co  from  the  class  $ computes  an  approximation  x^  such  that 

II  Al/<2 (^-0)11  S Ck||a^2H  ||  XjJI  c, 

»V“II  5 C*3/2l|Jkllc, 

llAvall s 5*ll  All  II  *JI c 

where  C = 3(5^  + 2)  + 0(C).  Furthermore  if  |j  Al/2  ||  ||  a\\  2=  ||  A1//2o||  then  the 
algorithm  cp  is  numerically  stable  in  the  A-norm  and  if  ||  A 2 1 1 ||  x, -a||  3: 

||  A (x^-a)  ||  then  the  algorithm  cp  is  stable  in  the  spectral  norm.  • 

Corollary  4.1  assures  that  the  algorithm  cp  computes  x^  with  the  relative 

3/2 

error  in  the  spectral  norm  depending  at  worst  on  (x  ' . The  residual  vector 
has  the  spectral  norm  of  order  at  most  £k.  We  repeat  that  the  algorithm  cp 
with  iterative  refinement  in  single  precision  is  numerically  stable  whenever 

3/2  2 

Cx  C < 1 and  well-behaved  whenever  CK  is  at  most  of  order  unity. 

We  now  give  an  example  of  an  algorithm  cd  which  satisfies  (4.6). 


4.5 


I Example  4 . 1 

■ - 

Let  and  be  the  computed  vectors  and  rk>  r^  the  corresponding 

residual  vectors.  Let  vk  - fl(Ark>  be  the  computed  vector  which  is  used  for 

the  computation  of  c,  . 

k 

We  propose  the  following  algorithm  for  the  computation  of  Let 

“1  • £l«V;k  - Vk»- 
“2  ■ - 'k  + Vk»- 

Thus  the  computation  of  and  does  not  require  further  matrix-vector  multi- 
plications. Repeating  a part  of  the  analysis  of  Section  3 it  is  possible  to 

I show  that 

”l  ’ <VA<V°»  + *V  16“!]  1 Cll  All  ||yk||  II^IICj 

“2  " * ®w2 ’ |6"2I  til  A||  ||?k||  ||^||c4 

where  — C^+l  and  — 2Cj+l.  From  this  we  get 
w 

' — - uk(l  -I-  5uk),  1*^1  f C||  All  ||  yk||  II^JI  (Cj/Jwjl  + c4/|w2|). 

This  suggests  the  following  algorithm  for  the  computation  of  u 

k’ 

^ 1E  til  *11  l|yk||  ||  (Cj/lwjl+c^lwjl)  < i 

\ - \ 

1 • 0 otherwise. 

Hence,  (4.6)  is  satisfied.  Note  that  u - 0 means  that  x,  ■ z » fi/C  _c  r 1 

k k+1  4k  ruxk  CkV 

is  obtained  by  one  step  of  the  steepest  descent  algorithm.  This  can  be  inter- 
preted as  the  initialization  of  the  eg  algorithm  from  the  vector  Xj^. 
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4.6 


It  may  also  be  observed  that  the  vectors  z^  and  y^  need  not  be  stored. 

One  step  of  the  algorithm  can  be  performed  having  five  vectors  x^,  x^  r^, 

^ and  vk  * Ar^  in  storage  and  using  two  matrix-vector  multiplications. 

We  have  performed  many  numerical  tests  using  this  algorithm.  In  most 
cases  the  algorithm  was  well-behaved  in  the  spectral  norm.  However,  in  a 
few  cases  (about  five  percent)  numerical  tests  experimentally  confirmed  the 
sharpness  of  the  error  bounds  in  Corollary  4.1.  ® 

2 

We  end  this  section  by  a remark  on  the  eg  algorithms  for  B ■ I and  B * A . 
Based  on  che  results  of  Section  3 and  assuming  that  the  computed  coefficient 

**k  * \(1  + 6uk)  where  l5ukl  s 1 and  uk " for  the 

computed  vectors  yk  and  z^,  it  is  possible  to  prove  that  there  exists  an 
index  k such  that  the  computed  x^  satisfies 

||  b1/2<v5||scMI  b1/2||  ll^llc. 

||  vS||*C||s1/2||  ||  b'1/2||  H^llc, 

II  «„-«»  s 6*11  a||  ||^||c 

for  a certain  constant  C “ C(n).  Note  that  for  B ■ I we  conclude  the  numerical 

stability  of  the  minimum  error  algorithm.  A detailed  analysis  for  the  minimal 

2 

residual  method,  B = A , may  be  found  in  Wieladek  [78]. 
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5.1 


5.  FINAL  COMMENTS 


We  have  shown  that  the  relative  error  of  the  computed  vector  by  a eg 

3/2 

algorithm  from  $ depends  at  worst  on  . Since  for  many  practical  cases 
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the  required  accuracy  is  larger  than  £x  ' this  is  a quite  satisfactory 
result. 

As  we  mentioned  before  we  have  not  succeeded  in  analyzing  classical  eg 
algorithms.  However,  we  believe  that  at  least  some  of  them  have  similar  num- 
erical properties. 

We  want  to  pose  another  problem  of  practical  interest  connected  with  the 
numerical  properties  of  eg  algorithms.  We  know  that  in  theory  the  sequence 
1x^3  approximates  the  solution  or  wi.th  the  best  possible  speed  of  convergence 
in  the  class  (2.5).  Is  this  still  true  in  the  presence  of  rounding  errors? 

It  is  important  to  know  the  speed  of  convergence  of  the  computed  sequence 
ix^}  and  to  see  how  much  of  the  theoretical  optimality  continues  to  hold  in  fl. 
We  observe  experimentally  that  Che  computed  sequence  initially  approximates  a 
at  least  as  fast  as  the  Chebyshev  iteration,  i.e., 

||  x^- Of | ^ s 2[  (a/k  - 1)/(«A  + l)]k||  xQ-a||  . Furthermore  in  many  cases  the  error 
||  x^-a|j^  is  significantly  less  than  the  above  bound.  Therefore  we  propose  the 
following  conjecture. 


for  all  k where  C ■ C(n). 
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Thus,  (5.1)  means  that  as  long  as  ||  is  greater  than  C*|| 
we  have  at  least  Chebyshev  speed  of  convergence.  For  large  k,  (5.1) 
numerical  stability  in  the  A-norm. 


3lkC* 

means 


A.  1 


APPENDIX 


We  describe  numerical  tests  of  the  Concus,  Golub  and  O'Leary  [76]  conjugate 
gradient  algorithm  defined  as  follows: 


Xq  - a given  approximation 


rk  " Axk  • b’ 


_ <vrk>  , 
k <vA?k> 


<rk*rk> 


-i 


k+1 


^k-l'V^  Ck-1 


w * 1 
1 1 


*k+l  " Xk-1  + Wk+l(ckrk  + \ ' Vl^  X-1  = 


for  k * 0,1 

— aT1  — * 

We  tested  this  algorithm  for  the  matrix  A = (I  * 2w  w )D(I  - 2w  w ) where 

w»  II  wll  * was  a vector  produced  by  a subroutine  which  generates  random 

numbers  and  D * diag  (X^,  , . . . , X.^)  was  a diagonal  matrix  with  Xj^  > 0.  We 

chose  n S [50,200]  for  different  distributions  of  varying  the  condition 
2 8 

number  from  10  to  10  . We  defined  X ^ as 

(i)  X^  ■ a + (1-a) (i-l)/(n- 1)  for  a positive  small  a, 

(ii)  Xj  ■ q,  X^  ■ a + (1-a) ( i— 2)/ (n-2) , i ■ 2, . . . ,n  where 
0 < q « a * 1, 

(iii)  X^  ■ q11  * for  some  q ^ 1. 


We  computed  j iterative  steps  until  the  limiting  accuracy  was  achieved. 
For  ill-conditioned  problems  with  n * 50,  j was  several  thousands.  As  we 


mentioned  in  the  Introduction  the  best  possible  computed  approximation  had 

3/2  - 14 

the  relative  error  of  order  £ x (where  £ equals  10  for  the  CDC  3600 

computer  used  in  the  experiments)  and  the  residual  vector  had  norm  of  order 

£k||  a||  . 

We  also  tested  the  algorithm  cp  described  in  Example  4.1  for  the  above 

examples.  In  most  cases  the  algorithm  cp  produced  the  residual  vectors  of  order 

£||  a||  . Thus  it  behaved  better  than  the  Concus,  Golub,  O'Leary  algorithm.  For 

a few  cases  with  eigenvalues  distributed  as  in  (iii)  and  with  initial  error 
n 

x„-o  ■ , c.C.  such  that  |c, | » |c.|  » ...  » |c  |,  the  algorithm  cp  produced 

0 4jJJ  '1'  1 2 1 n ' 

j-1 


the  residual  vectors  of  order  £x||  a||  . 
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